﻿//
///*
//#include "math.h"
//#include "stdio.h"
//#include "float.h"
//
//#include "ns.h"
//
//int holds(int n, hpdouble a)
//{
///*
//24.05.96
//Решение системы  A * X = B
//симм. положит. опр. матрица - метод кв.корня
//
//a(1:n*(n+1)/2) - верхняя треугольная часть
//возвращает :
//0 - o'k
//1 - нарушена пол.определенность
//*/
// double s;
// int k, ki, ik, jk;
// int nn, i, j, ij, i1;
// int ier;
//
// ier = 0;
///*..разложение A = L U */
//  nn = n + n;
//  ij = 1;
//  for( i=1; i<=n; i++ ) {
//    i1 = i - 1;
//    for( j=i; j<=n; j++ ) {
//      s = a[ij-1];
//      for( k=1; k<=i1; k++ ) {
//        ki = (k-1)*(nn-k)/2;
//        ik = ki + i;
//        jk = ki + j;
//        s -= a[jk-1]*a[ik-1];
//      }
//      if(j == i) {
///*..выход, если нарушено A <= 0 */
//        if(s <= MEPS) {
//          ier = 1;
//          goto lkex;
//        }
//        a[ij-1] = sqrt(s);
//      } else {
//        a[ij-1] = s/a[(i-1)*(nn-i)/2+i-1];
//      }
//      ij++;
//    }
// }
//lkex:
// return ier;
//}
//
//void hols(int n, hpdouble a, hpdouble b)
//{
///*
//24.05.96
//Решение системы  A * X = B
//симм. положит. опр. матрица метод кв.корня
//ч.2 - прямая и обратная подстановки
//
//a(0:n*(n+1)/2-1) - L'L разложение
//b(0:n-1)         - правая часть
//*/
//      double s;
//      int nn, j, k, jk, j1, nj, nk, k1, kj;
//
///* c..решение системы  L Y = B */
//      nn = n + n;
//      for( j=1; j<=n; j++ ) {
//      s = b[j-1];
//      j1 = j-1;
//      for( k=1; k<=j1; k++ ) {
//         kj=(k-1)*(nn-k)/2+j;
//         s -= a[kj-1]*b[k-1];
//      }
//         b[j-1] = s / a[(j-1)*(nn-j)/2+j-1];
//      }
///* c..решение системы  U X = Y      */
//      for( j=1; j<=n; j++ ) {
//         nj = n-j+1;
//         k1 = nj+1;
//         s=b[nj-1];
//         nk=n;
//         for( k=k1; k<=n; k++ ) {
//           jk = (nj-1)*(nn-nj)/2+nk;
//           s -= a[jk-1]*b[nk-1];
//           nk--;
//         }
//         b[nj-1] = s / a[(nj-1)*(nn-nj)/2+nj-1];
//      }
//      return;
//}
//
//void hols1(int n, hpdouble a, hpdouble b)
//{
///*
// 24.05.96
// Расчет новой строки L'L-разложения при пересчете матрицы Грама ( NSINEM )
//     a(0:n*(n+1)/2-1) - L'L разложение
//     b(0:n-1)         - новая строка
//*/
//      double s;
//      int nn, j, k, j1, kj;
///* c..решение системы  L Y = B */
//      nn = n + n;
//      for( j=1; j<=n; j++ ) {
//        s = b[j-1];
//              j1 = j-1;
//              for( k=1; k<=j1; k++ ) {
//                kj=(k-1)*(nn-k)/2+j;
//                s -= a[kj-1]*b[k-1];
//              }
//        b[j-1] = s / a[(j-1)*(nn-j)/2+j-1];
//      }
//      return;
//}
///**/
//
//double decomp(int n, hpdouble a, hpint ipvt, hpdouble work)
//{
//#define A(I_,J_) (a+(int)(I_)*(n)+(J_))
//
//int i, j, k, kb, km1, kp1, m, nm1;
//double anorm, ek, t, ynorm, znorm;
//double cond;
//
//ipvt[n-1] = 1;
//if(n == 1) goto L_80;
//nm1 = n - 1;
//anorm = 0.;
//for(j=1; j <= n; j++){
//  t = 0.;
//  for(i=1; i <= n; i++) t += fabs( *A(j-1,i-1) );
//  if( t > anorm ) anorm = t;
//}
//for(k=1; k <= nm1; k++){
//        kp1 = k + 1;
//        m = k;
//        for(i=kp1; i <= n; i++)
//           if( fabs( *A(k-1,i-1)) > fabs( *A(k-1,m-1)) ) m = i;
//        ipvt[k-1] = m;
//        if( m != k ) ipvt[n-1] =  -ipvt[n-1];
//        t =  *A(k-1,m-1);
//        *A(k-1,m-1) =  *A(k-1,k-1);
//        *A(k-1,k-1) = t;
//        if(t == 0.) goto L_35;
//        for(i=kp1; i <= n; i++) *A(k-1,i-1) =  - *A(k-1,i-1)/t;
//        for(j=kp1; j <= n; j++){
//                t =  *A(j-1,m-1);
//                *A(j-1,m-1) =  *A(j-1,k-1);
//                *A(j-1,k-1) = t;
//                if(t == 0.) goto L_30;
//                for(i=kp1; i <= n; i++)
//                         *A(j-1,i-1) =  *A(j-1,i-1) +  *A(k-1,i-1)*t;
//                L_30:
//                        ;
//        }
//        L_35:
//                ;
//}
//for(k=1; k <= n; k++){
//        t = 0.;
//        if( k == 1 ) goto L_45;
//          km1 = k - 1;
//          for(i=1; i <= km1; i++) t += *A(k-1,i-1)*work[i-1];
//        L_45:
//        ek = 1.;
//        if(t < 0.) ek = -1.;
//        if( *A(k-1,k-1) == 0. ) goto L_90;
//        work[k-1] =  - (ek+t)/ *A(k-1,k-1);
//                ;
//}
//for(kb=1; kb <= nm1; kb++){
//        k = n - kb;
//        t = 0.;
//        kp1 = k + 1;
//        for(i=kp1; i <= n; i++) t += *A(k-1,i-1)*work[k-1];
//        work[k-1] = t;
//        m = ipvt[k-1];
//        if( m == k ) goto L_60;
//        t = work[m-1];
//        work[m-1] = work[k-1];
//        work[k-1] = t;
//        L_60:
//                ;
//}
//ynorm = 0.;
//for(i=1; i <= n; i++) ynorm += fabs(work[i-1]);
//solve(n,a,work,ipvt);
//znorm = 0.;
//for(i=1; i <= n; i++) znorm += fabs(work[i-1]);
//cond = anorm*znorm/ynorm;
//if(cond < 1.) cond = 1.;
//return cond;
//L_80 :
//    cond = 1.;
//    if( *A(1-1,1-1) != 0. ) return cond;
//L_90 :
//    cond = 1.0e32;
//return cond;
//#undef A
//}
//
//void solve(int n, hpdouble a, hpdouble b, hpint ipvt)
//{
//#define A(I_,J_) (a+(int)(I_)*(n)+(J_))
//int i, k, kb, km1, kp1, m, nm1;
//double t;
//
//if( n == 1 ) goto L_50;
//nm1 = n - 1;
//for(k=1; k <= nm1; k++){
//        kp1 = k + 1;
//        m = ipvt[k-1];
//        t = b[m-1];
//        b[m-1] = b[k-1];
//        b[k-1] = t;
//        for(i=kp1; i <= n; i++) b[i-1] += *A(k-1,i-1)*t;
//}
//for(kb=1; kb <= nm1; kb++){
//        km1 = n - kb;
//        k = km1 + 1;
//        b[k-1] = b[k-1]/ *A(k-1,k-1);
//        t =  - b[k-1];
//        for(i=1; i <= km1; i++) b[i-1] += *A(k-1,i-1)*t;
//}
//L_50:
//  b[1-1] /= *A(1-1,1-1);
//return;
//#undef A
//}
//
///* nsinem 24.05.96 */
///**/
//void nsprm(int itr, int m1, int m2, int mk,
//           hpint ind, double alk,
//           hpdouble grm, hpdouble u, hpdouble v,
//           int ier, FILE *fl)
//{
//double elgrm;
//int i, ii, j, jj, m, mm, it;
//int ij;
//double s, si, sj, xx, yy;
//
//if(fl == NULL) return;
//
//if(itr == 0)
//  fprintf(fl,"NSINEM - NSPRM :\n\n Итер. Н.проек.   Н.прибл.   Шаг       Акт.  Номера :\n");
//
//m = m1 + m2;
//mm = m + m;
//
//yy = 0.;
//
//for(i=1; i <= mk; i++){
//        ii = ind[i-1];
//        if(ii > m) {
//          ii -= m2;
//          si = -1.;
//        } else si = 1.;
//        s = 0.;
//        for(j=1; j <= mk; j++){
//                jj = ind[j-1];
//                if( jj > m ) {
//                  jj -= m2;
//                  sj = -1.;
//                } else sj = 1.;
//    if(ii<jj) ij = (int)(mm-ii)*(ii-1)/2 + jj;
//    else ij = (int)(mm-jj)*(jj-1)/2 + ii;
//                elgrm = grm[ij-1];
//                s += si*sj*elgrm*v[j-1];
//        }
//        yy += s*v[i-1];
//}
//
//if(yy > 0.) yy = sqrt(yy);
//
//xx = 0.;
//for(i=1; i <= m; i++){
//        s = 0.;
//        for(j=1; j <= m; j++){
//    if(i<j) ij = (int)(mm-i)*(i-1)/2 + j;
//    else ij = (int)(mm-j)*(j-1)/2 + i;
//                elgrm = grm[ij-1];
//                if(j <= m1) s += elgrm*u[j-1];
//                else s += elgrm*(u[j-1]-u[j+m2-1]);
//        }
//        if(i <= m1) xx += s*u[i-1];
//        else xx += s*(u[i-1]-u[i+m2-1]);
//}
//
//if(xx > 0.) xx = sqrt(xx);
//
//it = (itr > 0)? itr : -itr;
//
//fprintf(fl,"\n%4d %10.3le %10.3le %10.3le %4d  ", it, yy, xx, alk, mk );
//
//ii=0;
//for(i=1;i<=mk;i++) {
//  if(ii >= 10) {
//    ii = 0;
//    fprintf(fl,"\n %44s"," " );
//  }
//  fprintf(fl,"%3d ", ind[i-1] );
//  ii++;
//}
//
//if(ier == 1) fprintf(fl,"\n\n ПРЕРЫВАНИЕ HOLDS\n\n");
//if(ier == 3) fprintf(fl,
//"\n\n ОШИБКА В ВХОДНЫХ ДАННЫХ ЛИБО ПЛОХАЯ ОБУСЛОВЛЕННОСТЬ МАТРИЦЫ ГРАМА\n\n");
//if(itr < 0) {
//    if(ier == 0) fprintf(fl,"\n\n ОПТИМАЛЬНЫЕ МНОЖИТЕЛИ :\n\n");
//    if(ier == 2) fprintf(fl,"\n\n ЛИМИТ ИТЕРАЦИЙ\n\n");
//    ii=0;
//    for(i=1; i<=mk; i++) {
//       if(ii >= 5) {
//                ii = 0;
//                fprintf(fl,"\n");
//       };
//       fprintf(fl,"%10.4lg ", v[i-1] );
//       ii++;
//    };
//}
//
//return;
//}
//
//void compr(int mkp, int idx, hpdouble grm)
//{
//  int i, j, jw, jc, mk, mkp2, mk2;
//  int ijp, ij;
//
//  mk = mkp - 1;
//  mkp2 = 2*mkp;
//  mk2 = 2*mk;
//
//  for(i=1;i<=mk;i++) {
//    jw = i;
//    jc = i;
//    if(i > idx) jc++;
//    for(j=jc;j<=mkp;j++)
//      if(j != idx) {
//        if(i<j) ijp = (int)(mkp2-i)*(i-1)/2 + j;
//        else    ijp = (int)(mkp2-j)*(j-1)/2 + i;
//        if(i<jw) ij = (int)(mk2-i)*(i-1)/2 + jw;
//        else     ij = (int)(mk2-jw)*(jw-1)/2 + i;
//        grm[ij-1] = grm[ijp-1];
//        jw++;
//      };
//  };
//
//  return;
//}
//
//void prn_m(int mk, hpdouble grm, FILE *fl)
//{
//  int i, j, mk2;
//  int ij;
//  double z;
//
//  z = 0.;
//  mk2 = 2*mk;
//
//  fprintf(fl,"\n\n");
//  for(i=1;i<=mk;i++) {
//    fprintf(fl,"\n");
//    for(j=1;j<=mk;j++) {
//        if(i<=j) {
//          ij = (int)(mk2-i)*(i-1)/2 + j;
//          fprintf(fl,"%10.3le ",grm[ij-1]);
//        } else {
//          ij = (int)(mk2-j)*(j-1)/2 + i;
//          fprintf(fl,"%10.3le ",z);
//        };
//    };
//  };
//
//  return;
//}
//
//
//double rsqr(double x, double y)
//{
// double u, v, w;
//
// x = fabs(x);  y = fabs(y);
// v = (x >= y)?x:y;
// u = (x <= y)?x:y;
// if(v == 0.) w = 0.;
// else {
//   w = u/v;
//   w = v*sqrt(1.+w*w);
// };
// return w;
//}
//
//int nsinem(int m1, int m2, int *mkp, hpint ind,
//           hpdouble grm, hpdouble grk,
//           hpdouble v, hpdouble u, hpdouble c,
//           int *litp, FILE *fl)
//{
///*
//   (C) ГOPБУHOB В.К. MAЙ 87 - ДEK 88, КОХАНОВСКИЙ И.И. МАЙ 96
//   HOPMAЛЬHOE PEШEHИE CИCTEMЫ
//       A(I)*X = C(I), 1<=I<=M1
//       MOD(A(M1+I)*X-Y(I)) <= DEL(I), 1<=I<=M2
//       HAЧAЛЬHAЯ TOЧKA - HOPM. PEШ.  A*X=C
//        M=M1+M2 (PAЗMEPHOCTЬ ЗAДAЧИ)
//        MK-AKTИBHЫX OГPAHИЧEHИЙ(<=M)
//        NK-HEAKTИBHЫX (2*M2)
//        IND(1:MK)-HOMEPA AKT. OГPAH.
//        IND(M+1:M+NK)-HOMEPA HEAKTИBHЫX             ! 3M
//        GRK-MATPИЦA ГPAMA AKTИBHЫX OГPAHИЧEHИЙ
//        N=M+M2
//        V(1:MK)-KOЭФФ. PAЗЛOЖEHИЯ X=V*A             ! 2M
//        V(M+1:2M)- ИСПОЛЬЗУЕТСЯ ПРИ пересчете матрицы Грама
//        U(1:N)-KOЭФФ. PAЗЛOЖ. TOЧKИ XK              ! 4M
//        U(N+1,2*N)-KOЭФФ. PAЗЛOЖ HAПPABЛEHИЯ YK-XK
//        C(M1+I)=Y(I)+DEL(I),1<=I<=M2                ! 2M
//        C(M+I)=-Y(I)+DEL(I),1<=I<=M2
//        LIT-ЛИMИT ИTEPAЦИЙ
//        LPR-ЛИMИT ПEЧATИ
//     ПOДПPOГPAMMЫ:
//         HOLDS-PEШEHИE CЛAУ
//         NSPRM-ПPOMEЖУTOЧHAЯ ПEЧATЬ
//
// ier = 0 - O'K
//       1 - авар.выход HOLDS - вырожденность матрицы Грама
//       2 - лимит итераций
//       3 - ошибка в входных данных,
//           либо плохая обусловленность матрицы Грама
//*/
//const double ALKMAX=1.e10;
//
//double eps, seps;
//double elgrm, eldiag;
//int itr, ier, lit;
//int i, j, m, n, mk;
//int ii, jj, im, imk, m11, ni, nj, nk, nm2, imov, mm, i_del=0;
//int ij, ijk;
//double ali, alk, di, ei, one, s, si, sj;
//int mko, mko2, mk2;
//int F_ADD=FALSE, F_DEL=FALSE, F_OPT=FALSE, fe;
//
//  ier = 0;
//  lit = *litp;
//  eps=MEPS;
//  seps = sqrt(eps);
//  one = 1. - seps;
///**/
//  m = m1 + m2;
//  mm = m + m;
//  mk = m1;
//  m11 = m1 + 1;
//  n = m + m2;
//  nm2 = n + m2;
//  nk = n - mk;
//
//  for(j=1; j <= nk; j++) ind[j+m-1] = j + m1;
//
//  for(i=1; i <= m; i++){
//    u[i-1] = v[i-1];
//    u[i+m-1] = 0.;
//    if( i > m1 ) ind[i-1] = 0;
//    v[i-1] = 0.;
//  };
//  if(fl != NULL) nsprm(0,m1,m2,mk,ind,0.,grm,u,v,ier,fl);
//
///*                                      ИTEPAЦИИ AЛГOPИTMA               */
//  for(itr=1; itr <= lit; itr++){
//        if( mk == 0 ) goto L_MK0;
///*                            ФOPMИPOBAHИE ГPAHИ     */
//
///* здесь mk >= 1 */
//
//        if(!F_ADD && !F_DEL) {
//          ijk = 0;
//          for(j=1; j <= mk; j++){
//                jj = ind[j-1];
//                v[j-1] = c[jj-1];
//
//                if(jj > m) {
//                  jj -= m2;
//                  sj = -1.;
//                } else sj = 1.;
//                for(i=j; i <= mk; i++ ){
//                   ii = ind[i-1];
//
//                   if(ii > m)  {
//                     ii -= m2;
//                     si = -1.;
//                   } else si = 1.;
//       if(ii<jj) ij = (int)(mm-ii)*(ii-1)/2 + jj;
//       else  ij = (int)(mm-jj)*(jj-1)/2 + ii;
//                   ijk++;
//                   grk[ijk-1] = si*sj*grm[ij-1];
//                };
//          };
//
//          ier = holds( mk, grk );
//          if(ier > 0) {
//      ier = -ier;
//      goto LK_EX;
//    }
//        };
//
//  if(F_DEL) {
//    int mkp;
//
//    F_DEL = FALSE;
//    mkp = mk + 1;
//
//          for(j=1; j<=mk; j++){
//                  jj = ind[j-1];
//                  v[j-1] = c[jj-1];
//    };
//
//    if(i_del == mkp) compr(mkp, i_del, grk);
//    else {
//      int t, k, mkp2, iw;
//      int  ij1, ij2;
//      double c, s, g1, g2, rg;
//
//      mkp2 = 2*mkp;
//      for(t=0; t<=mkp-1-i_del; t++) {
//        iw = i_del + t;
//        ii = iw + 1;
//        jj = iw;
//        ij1 = (int)(mkp2-jj)*(jj-1)/2 + ii;  /* ii >= jj */
//        g1 = grk[ij1-1];
//
//        jj++;
//        ij2 = (int)(mkp2-jj)*(jj-1)/2 + ii;  /* ii >= jj */
//        g2 = grk[ij2-1];
//
//        rg = rsqr(g1,g2);
//        if( rg == 0. ) {
//          ier = -1;
//          goto LK_EX;
//        };
//        c = g1/rg;
//        s = g2/rg;
//
//        for(k=1; k<=mkp-i_del-t; k++) {
//          ii = iw + k;
//          jj = iw;
//          ij1 = (int)(mkp2-jj)*(jj-1)/2 + ii;   /* ii >= jj */
//          g1 = grk[ij1-1];
//
//          jj++;
//          ij2 = (int)(mkp2-jj)*(jj-1)/2 + ii;   /* ii >= jj */
//          g2 = grk[ij2-1];
//
//          grk[ij1-1] =  c*g1 + s*g2;
//          grk[ij2-1] = -s*g1 + c*g2;
//        };
//      };
//      compr(mkp, i_del, grk);
//    };
//  };   /* FDEL */
//
//  if(F_ADD) {
///* добавлено ограничение, F_ADD=TRUE */
///* номер добавленного ограничения - последний в массиве IND */
//          F_ADD = FALSE;
//          mko = mk - 1;
//          mko2 = 2 * mko;
//          ijk = 0;
//          for(j=1; j <= mk; j++){
//                 jj = ind[j-1];
//                 v[j-1] = c[jj-1];
//
//                 if(jj > m) {
//                   jj -= m2;
//                   sj = -1.;
//                 } else sj = 1.;
///* grm(ik,j) j-из активного набора; ik-новое ограничение */
//                 ii = ind[mk-1];
//                 if(ii > m)  {
//                   ii -= m2;
//                   si = -1.;
//                 } else si = 1.;
//     if(ii<jj) ij = (int)(mm-ii)*(ii-1)/2 + jj;
//     else  ij = (int)(mm-jj)*(jj-1)/2 + ii;
//                 ijk++;
//                 v[m+ijk-1] = si*sj*grm[ij-1];
//          };
//          hols1( mk-1, grk, &(v[m]) );
///* расчет нового диагонального элемента */
//          s=0.;
//          for(i=0; i<mk-1; i++) s += v[m+i]*v[m+i];
//          s = v[m+mk-1] - s;
//          if(s <= eps) {
///* малый диагональный элемент, текущая матрица Грама вырождена */
//            ier = -1;
//            goto LK_EX;
//          } else {
//            eldiag = sqrt(s);
///* раздвигаем матрицу Грама..             */
//            for(i=1; i<=mko; i++) {
//              ni = mko - i + 1;
//              imov = ni - 1;
//              for(j=ni; j<=mko; j++) {
//                nj = mko - j + ni;
//                      ij = (mko2-ni)*(ni-1)/2 + nj;
//                      grk[ij+imov-1] = grk[ij-1];
//              };
//            };
///*.. и дописываем в нее новый столбец */
//            mk2 = 2*mk;
//            for(i=1; i<=mk; i++) {
//              ij = (mk2-i)*(i-1)/2 + mk;
//              if(i != mk) grk[ij-1] = v[m+i-1];
//              else grk[ij-1] = eldiag;
//            }
//          }
///* конец пересчета */
//        }
///*                            BЫЧИCЛEHИE ПPOEKЦИИ                        */
//        hols(mk, grk, v);
//
//        for(i=1; i <= mk; i++) {
//          ii = ind[i-1];
//          u[ii+n-1] = v[i-1] - u[ii-1];
//        }
///*                                  ПEPEXOД MK=0, ДOПУCTИMOCTЬ ПPOEKЦИИ  */
//        L_MK0:
//
//        for(j=1; j <= nk; j++){
//          jj = ind[j+m-1];
//          u[jj+n-1] =  -u[jj-1];
//        };
//
//        alk = ALKMAX;
///*                                                AHAЛИЗ OГPAHИЧEHИЙ     */
//        for(i=m11; i <= m; i++) {
//
//    fe = FALSE;
//                for(j=m11; j <= mk; j++)
//                   if(i == ind[j-1] || i == ind[j-1] - m2) {
//         fe = TRUE;
//         break;
//       };
//    if( fe ) continue;
//
//                di = 0.;
//                ei = 0.;
//                for(j=1; j <= m; j++) {
//      if(i<j) ij = (int)(mm-i)*(i-1)/2 + j;
//      else  ij = (int)(mm-j)*(j-1)/2 + i;
//        elgrm = grm[ij-1];
//                        if(j <= m1) {
//                          di += elgrm*u[j-1];
//                          ei += elgrm*u[n+j-1];
//                        } else {
//                          di += elgrm*(u[j-1]-u[m2+j-1]);
//                          ei += elgrm*(u[n+j-1]-u[nm2+j-1]);
//                        };
//                };
//
//                if(fabs(ei) <= eps) continue;
//
//                if(ei < -eps) {
//                  im = i + m2;
//                  ali = -(c[im-1]+di)/ei;
//                } else {
//                  im = i;
//                  ali = (c[im-1]-di)/ei;
//                };
//
//                if(alk > ali) {
//                  alk = ali;
//                  imk = im;
//                };
//        };
//
//        if(fl != NULL) nsprm(itr,m1,m2,mk,ind,alk,grm,u,v,ier,fl);
//
//        if(alk < one) {
///*   ПPOEKЦИЯ HEДOПУCTИMAЯ  */
//          F_ADD = TRUE;
//          if(alk < -seps) {
//      ier = -3;
//      goto LK_EX;
//    }
//          if(alk < 0.) alk = 0.;
//          for(j=1; j <= n; j++) u[j-1] += alk*u[n+j-1];
///**/
//    if(itr == lit) goto L_LIT;
///**/
//
//          mk++;
//          ind[mk-1] = imk;
//          nj = nk;
//          for(j=1; j <= nj; j++) {
//                  if( ind[j+m-1] != imk )       continue;
//                  ind[j+m-1] = ind[nk+m-1];
//                  ind[nk+m-1] = 0;
//                  nk--;
//      break;
//          };
//  } else {
///*     ПPOEKЦИЯ ДOПУCTИMAЯ    */
///**/
//    if(itr == lit) goto L_LIT;
///**/
//          for(i=1; i <= mk; i++) u[ind[i-1]-1] = v[i-1];
//          if(mk == m1) {
//      F_OPT = TRUE;
//      goto L_OPT;
//    };
//
//          for(j=1; j <= nk; j++) u[ind[j+m-1]-1] = 0.;
//
//    F_OPT = TRUE;
//          ni = mk;
//          for(i=m11; i <= ni; i++) {
//                  if(v[i-1] <= 0.) continue;
//      F_OPT = FALSE;
//      F_DEL = TRUE;
//      i_del = i;
//                  nk++;
//                  ind[nk+m-1] = ind[i-1];
//      for(j=i; j<=mk-1; j++) ind[j-1] = ind[j];
//                  ind[mk-1] = 0;
//                  mk--;
//                  break;
//        };
//  };
//  if(F_OPT) goto L_OPT;
//  };
///*     OБPAБOTKA PEЗУЛЬTATA   */
//L_LIT :         ;
///*     ЛИМИТ ИТЕРАЦИЙ          */
//  ier = 2;
//L_OPT :         ;
///*     TOЧKA OПTИMAЛЬHAЯ      */
//  itr = -itr;
//  if(fl != NULL) nsprm(itr,m1,m2,mk,ind,alk,grm,u,v,ier,fl);
//
//  itr = (itr>0) ? itr : -itr;
//
//  if(F_OPT) {
//    for(j=1; j <= n; j++) u[j-1] = 0.;
//    for(i=1; i <= mk; i++)
//      if(ind[i-1] > m) {
//        ind[i-1] -= m2;
//        u[ind[i-1]-1] = -v[i-1];
//      } else {
//        u[ind[i-1]-1] = v[i-1];
//      };
//    ind[mk] = 0;
//    for(i=1; i<=m; i++) v[i-1] = u[i-1];
//  } else {
///* ier = 2 */
//    for(i=m11; i <= m; i++) u[i-1] -= u[m2+i-1];
//    ind[mk] = 0;
//    for(i=1; i<=m; i++) v[i-1] = u[i-1];
//  };
//
//LK_EX :
///*     АВАРИЙНОЕ ЗАВЕРШЕНИЕ      */
//        if(ier < 0) {
//    ier = -ier;
//    if(fl != NULL) nsprm(itr,m1,m2,mk,ind,alk,grm,u,v,ier,fl);
//  };
//
//  *mkp = mk;
//  *litp = itr;
//
//  return ier;
//}
///**/
///* evd */
///**/
//
//int tred2( int n, hpdouble a, hpdouble d, hpdouble e,
//                            hpdouble z, double tol );
//int imtql2( int n, hpdouble z, hpdouble d, hpdouble e, double eps );
//
//int evdp( int ncall, int n, hpdouble pD,
//              double tola, double tolb, double cond, int irang )
///*
// ier = 0 - o"k
//       -1 - ошибка при записи в файл
//*/
//{
//int ier=0;
//string fn="nsrep_ev.c";
//int i, j;
//FILE* fl;
//
//if( ncall == 0 ) {
//   if( (fl = fopen( fn, "wt" )) == NULL ) {
//    ier = -1;
//    goto lkex;
//   }
//   fprintf(fl,"EVD - EVDP :\n\n");
//   fprintf(fl,"Размерность     : %4d\n",n);
//   fprintf(fl,"\nСобственные числа :\n\n");
//   j = 0;
//   for( i=0; i<n; i++) {
//      fprintf(fl," %3d %10.4lg",i+1,pD[i]);
//      j++;
//      if( j>4 ) { j=0; fprintf(fl,"\n"); }
//   }
//} else {
//   if( (fl = fopen( fn, "at" )) == NULL ) {
//    ier = -1;
//    goto lkex;
//   }
//   fprintf(fl,"\ntola = %10.4lg\n",tola);
//   fprintf(fl,"tolb = %10.4lg\n",tolb);
//   fprintf(fl,"irang = %3d\n",irang);
//   fprintf(fl,"cond  = %10.4lg\n",cond);
//}
//
//lkex :
//
//fclose(fl);
//
//return ier;
//}
//
//int evdt( int n, hpdouble ptola )
///* Формирование (или считывание) оценки tola*/
//{
//int ier=0;
//string fn="evd_t.c";
//double tola, eps;
//FILE* fl;
//char slpr[81];
//
//  if( (fl = fopen( fn, "rt" )) != NULL ) {
//    if( fgets( slpr, 81, fl ) != NULL ) fscanf( fl, "%le", &tola );
//    if( tola < 0. || tola > 1. ) tola = 0.;
//    fclose(fl);
//  } else tola = 0.;
//
//  if( tola == 0. ) {
///* оценка по умолчанию  : || aij - aij~ || <= n**2*||A||*2**(1-t) */
//   eps = pow(2,1-DBL_MANT_DIG); 
//   tola = n*eps;
//  }
//*ptola = tola;
//return ier;
//}
//
//int evdd( int n, hpdouble a, hpdouble d )
//{
///* .. 21.10.93 .. 08.12.93 "C" версия 2D_APP */
///*     Вычисление спектрального разложения симметричной матрицы A : */
///*     A = U * D * U'        (*) */
///*     D - диагональная матрица, элементы которой - собственные числа */
///*         матрицы A, упорядоченные по убыванию; */
///*     U - ортогональная матрица, столбцы которой - соответственно */
///*         собственные вектора A. */
///*     Параметры : */
///*     n        - размерность задачи; */
///*     a(n,n) - на входе : */
///*                исходная симметричная матрица A; */
///*                на выходе : */
///*                ортогональная матрица U из (*); */
///*     d(2*n)   - на выходе : */
///*                d(1:n) - собственные числа матрицы A (диагональ D); */
///*                d(n+1:2*n) - рабочий массив ( наддиагональ трех- */
///*                диагональной матрицы ортогонально эквивалентной A ); */
///*     ier      - на выходе : */
///*                0 - O'K; */
///*                1 - неудача. Для определения собственного числа */
///*                    требуется > limitr итераций QR алгоритма. */
///*.. */
//
//double eps=MEPS;
//int ier=0;
//
///*..Приведение матрицы A к трехдиагональной (отражения Хаусхолдера) */
//ier=tred2(n,a,d,(d+n),a,eps*eps);
///*..Определение собственных чисел и векторов A (QL алгоритм с неявным сдвигом) */
//if( ier==0 ) ier = imtql2(n,a,d,(d+n),eps);
//
//return ier;
//}
//
//int evds( int n, hpdouble a, hpdouble d, hpdouble b, hpdouble x,
//          double tola, double tolb, hpint irang )
//{
///* ****** 2-а отличия от ФОРТРАН версии */
///*.. 06.01.94 */
///*     Решение системы линейных уравнений на основе спектрального разложения. */
///*     Матрица системы A - симметричная, положительно определенная. */
///*     A * X = B,            (*) */
///*     A = U * D * U'        (**) */
///*     D - диагональная матрица, элементы которой - собственные числа */
///*         матрицы A, упорядоченные по убыванию; */
///*     U - ортогональная матрица, столбцы которой - соответственно */
///*         собственные вектора A. */
///*     т.о. */
///*       D * (U' * X) = U' * B, */
///*       Y = U' * X, */
///*       D * Y = U' * b, */
///*       X = U * Y. */
///*     Параметры : */
///*     n        - размерность задачи; */
///*     a(n,n) - на входе : */
///*                ортогональная матрица U из (**); */
///*     d(2*n)     - на входе : (в отличии от ФОРТРАН версии) */
///*                d(0:n-1) - собственные числа матрицы A (диагональ D), */
///*                d(n:2*n-1) - рабочий массив */
///*     b(n)     - на входе : */
///*                правая часть системы ограничений; */
///*     x(n)     - на выходе : */
///*                решение системы (*); */
///*     tola     - на входе : */
///*                норма матрицы погрещностей : */
///*                собственные числа < tola считаются нулевыми */
///*     tolb     - оценка ср.квадратичной погрешности правой части */
///*     irang    - r-ранг A    */
///*     ier      - на выходе : */
///*                0 - O'K; */
///*                2 - матрица системы не положительноопределенная */
///*.. */
//#define A(I_,J_) (a+(int)(I_)*(n)+(J_))
//int i, j, k;
//int kex;
//double eps1, eps2, s, tol, tolb2;
//double eps=MEPS;
//int f_neg=0;
//int ier = 0;
//
//tola = fabs(tola);
//eps1 = eps*sqrt( (double)n );
//eps2 = sqrt(eps1);
//for(j=1; j <= n; j++ ) {
//   if( tola > 0. ) {
//     if( fabs(d[j-1]) <= eps1 ) d[j-1] = 0.;
//     if( d[j-1] < 0. && d[j-1] >= -eps2 ) d[j-1] = 0.;
//   }
//   if( d[j-1] < -eps2 ) f_neg = 1;
//}
///*.. */
//if( f_neg == 1 ){
//        /*.. матрица системы не положительноопределенная */
//  ier = 2;
//  goto L_99;
//}
///*.. */
///*.. B~ = U' * B */
//for(j=1; j <= n; j++ ){
//   s = 0.;
//   for(i=1; i <= n; i++ ) s += *A(i-1,j-1)*b[i-1];
//   x[j-1] = s;
//}
//
///*.. Сглаживание возмущений правой части */
//if( tolb > 0. ){
//  tolb2 = tolb*tolb;
//  kex = n + 1;
//  s = 0.;
//  for(k=n; k >= 1; k--  ) {
//      s += x[k-1]*x[k-1];
//      if( s > tolb2 ) {
//        kex = k + 1;
//        break;
//      }
//  }
//  if( kex < n ) for(k=kex; k <= n; k++ ) d[k-1] = 0.;
//}
///*.. D * Y = B~ */
//
//tol = (tola>tolb)?tola:tolb;
//*irang = 0;
//for(j=1; j <= n; j++ ){
//  if( d[j-1] > tol ) {
//    x[j-1] = x[j-1]/d[j-1];
//    *irang += 1;
//  } else {
//    x[j-1] = 0.;
///* ******  отличие N2 от ФОРТРАН версии */
//    d[j-1] = 0.;
//  }
//}
//
///*.. X = U * Y */
//
//for(i=1; i <= n; i++ ){
//   s = 0.;
//   for(j=1; j <= n; j++ ) s += *A(i-1,j-1)*x[j-1];
///* ******  отличие N1 от ФОРТРАН версии */
//   d[n+i-1] = s;
//}
///*.. */
//
//for(j=1; j <= n; j++ ) x[j-1] = d[n+j-1];
//
///*.. */
//
//L_99:
//
//return ier;
//
//#undef A
//}
//
//int tred2( int n, hpdouble a, hpdouble d, hpdouble e,
//                            hpdouble z, double tol )
//{
//#define A(I_,J_) (a+(int)(I_)*(n)+(J_))
//#define Z(I_,J_) (z+(int)(I_)*(n)+(J_))
///*     14.10.93 - */
///*     20.10.93 */
///* Приведение симметрической матрицы А к трехдиагональной форме */
///* с помощью преобразования Хаусхолдера с использованием лишь   */
///* диагональных и поддиагональных элементов. Исходная матрица А */
///* записана в массиве А[1:n,1:n]. Диагональные элементы результиру- */
///* ющей матрицы накапливаются в массиве D[1:n], а поддиагональные */
///* в последних (n-1) ячейках массива e[1:n] ( e[1]=0 ). Матрица */
///* преобразования заносится в массив Z[1:n,1:n]. Массив А сохраня- */
///* ется, если фактические значения формальных параметров A и Z не */
///* совпадают. */
///* Уилкинсон Дж.Х., Райнш С. Справочник алгоритмов на языке АЛГОЛ. */
///* Линейная алгебра. ( c.194 ). */
//
//int i, j, k, l;
//double f, g, h, hh;
//int ier=0;
//
//for(i=1; i <= n; i++ ) {
//        for(j=1; j <= i; j++ ) *Z(i-1,j-1) =  *A(i-1,j-1);
//}
//
//for(i=n; i >= 2; i-- ){
//        l = i - 2;
//        f =  *Z(i-1,i-1-1);
//        g = 0.;
//        for(k=1; k <= l; k++ )  g += *Z(i-1,k-1)* *Z(i-1,k-1);
//        h = g + f*f;
//        /*..переобразование Хаусхолдера не выполняется и происходит выход к метке 999 */
//        /*..если параметр g слишком мал чтобы гарантировать ортогональность */
//        if( g <= tol ) {
//            e[i-1] = f;
//            h = 0.;
//            goto L_999;
//        }
//        l++;
//
//        if( f >= 0. ) e[i-1] =  -sqrt(h);
//        else e[i-1] = sqrt(h);
//
//        g = e[i-1];
//        h -= f*g;
//        *Z(i-1,i-1-1) = f - g;
//        f = 0.;
//        for(j=1; j <= l; j++ ){
//                *Z(j-1,i-1) =  *Z(i-1,j-1)/h;
//                g = 0.;
//                /*.. формирование A*u */
//                for(k=1; k <= j; k++ )  g += *Z(j-1,k-1)* *Z(i-1,k-1);
//                for(k=j+1; k <= l; k++ ) g += *Z(k-1,j-1)* *Z(i-1,k-1);
//                /*.. формирование p */
//                e[j-1] = g/h;
//                f +=  g* *Z(j-1,i-1);
//        }
//        /*..формирование K */
//        hh = f/(h+h);
//        /*..преобразование A */
//        for(j=1; j <= l; j++ ){
//                f =  *Z(i-1,j-1);
//                e[j-1] -= hh*f;
//                g = e[j-1];
//                for(k=1; k <= j; k++ )
//                   *Z(j-1,k-1) =  *Z(j-1,k-1) - f*e[k-1] - g* *Z(i-1,k-1);
//
//        }
//L_999:
//                ;
//        d[i-1] = h;
//}
///*.. end do i.. */
//
//d[0] = e[0] = 0.;
//
///*..накопление матриц преобразований */
//
//for(i=1; i <= n; i++ ){
//        l = i - 1;
//        if( d[i-1] != 0. ) {
//                for(j=1; j <= l; j++ ){
//                        g = 0.;
//                        for(k=1; k <= l; k++ )
//                           g += *Z(i-1,k-1)* *Z(k-1,j-1);
//                        for(k=1; k <= l; k++ )
//                           *Z(k-1,j-1) -= g* *Z(k-1,i-1);
//                 }
//        }
//        /*.. */
//        d[i-1] =  *Z(i-1,i-1);
//         *Z(i-1,i-1) = 1.;
//        for(j=1; j <= l; j++ ){
//                 *Z(i-1,j-1) = 0.;
//                 *Z(j-1,i-1) =  *Z(i-1,j-1);
//        }
//}
//
//return ier;
//
//#undef Z
//#undef A
//}
//
//
//int imtql2( int n, hpdouble z, hpdouble d, hpdouble e, double eps )
//{
//#define Z(I_,J_) (z+(int)(I_)*(n)+(J_))
///*ccc */
///*     14.10.93 - */
///*     20.10.93   */
///*ccc */
//
//int i, ia, j, k, m;
//int its, ier=0;
//double c, h, p, q, s, t, u, wrk;
//double tmpZ;
//int limitr=300;
//
//for(i=2; i <= n; i++ ){
//        e[i-1-1] = e[i-1];
//}
//
//e[n-1] = 0.;
//k = n - 1;
//for(j=1; j <= n; j++ ){
//        its = 0;
//        /*.. выделение пренебрежимо малого по величине поддиагонального элемента */
//        /*..test : */
//L_10:
//                ;
//
//        for(m=j; m <= k; m++ ) {
//                /*.. пеpеход к cont1 */
//                if( fabs(e[m-1]) <= eps*(fabs(d[m-1])+fabs(d[m+1-1])) )
//                        goto L_20;
//        }
//        m = n;
//       /*.. cont1 : */
//L_20:
//                ;
//        u = d[j-1];
//        if( m != j ){
//                if( its == limitr ) {
//                        ier = 1;
//                        goto L_99;
//                };
//                its++;
//                /*..формирование сдвига */
//                q = (d[j+1-1]-u)/(2.*e[j-1]);
//                t = sqrt(1.+q*q);
//
//                if( q < 0. ) wrk = q - t;
//                else wrk = q + t;
//
//                q = d[m-1] - u + e[j-1]/wrk;
//                u = 0.;
//                c = 1.;
//                s = c;
//
//                for(i=m-1; i >= j; i-- ){
//                        p = s*e[i-1];
//                        h = c*e[i-1];
//                        if( fabs(p) >= fabs(q) ){
//                                c = q/p;
//                                t = sqrt(c*c+1.0);
//                                e[i+1-1] = p*t;
//                                s = 1.0/t;
//                                c = c*s;
//                        } else {
//                                s = p/q;
//                                t = sqrt(s*s+1.0);
//                                e[i+1-1] = q*t;
//                                c = 1.0/t;
//                                s = s*c;
//                        }
//                        q = d[i+1-1] - u;
//                        t = (d[i-1]-q)*s + 2.0*c*h;
//                        u = s*t;
//                        d[i+1-1] = q + u;
//                        q = c*t - h;
//                        /*.. Вычисление собственного вектора */
//						
//                        for(ia=1; ia <= n; ia++ ){
//                                tmpZ =  *Z(ia-1,i-1);
//                                p =  *Z(ia-1,i+1-1);
//                                *Z(ia-1,i+1-1) = s*tmpZ + c*p;
//                                *Z(ia-1,i-1) = c*tmpZ - s*p;
//                        }
//                }
//                /*..        enddo i */
//                d[j-1] -= u;
//                e[j-1] = q;
//                e[m-1] = 0.;
//
//                /*.. goto test.. */
//                goto L_10;
//
//                }
//        /*..      endif m.ne.j */
//}
///*..   enddo j */
//
///*.. размещение собственных значений и векторов в порядке убывания */
//
//for(i=1; i <= n; i++ ){
//        k = i;
//        p = d[i-1];
//        for(j=i+1; j <= n; j++ ) {
//                if( d[j-1] > p ) {
//                        k = j;
//                        p = d[j-1];
//                }
//        }
//        if( k != i ){
//                d[k-1] = d[i-1];
//                d[i-1] = p;
//                for(j=1; j <= n; j++ ){
//                        p =  *Z(j-1,i-1);
//                         *Z(j-1,i-1) =  *Z(j-1,k-1);
//                         *Z(j-1,k-1) = p;
//                }
//        }
//}
//
//L_99:
//
//return ier;
//
//#undef Z
//}
//*/
//
// 
